\begin{tiny}
\begin{verbatim}
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "gpu_galois.h"

/* Global variables for regions to
   avoid allocating separately on each 
   call */
unsigned char *c_ur1_d = NULL;
unsigned char *c_ur2_d = NULL;
unsigned short *s_ur1_d = NULL;
unsigned short *s_ur2_d = NULL;
unsigned int *l_ur1_d = NULL;
unsigned int *l_ur2_d = NULL;

/* Globals for primitive polynomials,
   one for device and one for host. */
int *prim_poly_d = NULL;

static int prim_poly[33] = 
{ 0, 
	/*  1 */     1, 
	/*  2 */    07,
	/*  3 */    013,
	/*  4 */    023,
	/*  5 */    045,
	/*  6 */    0103,
	/*  7 */    0211,
	/*  8 */    0435,
	/*  9 */    01021,
	/* 10 */    02011,
	/* 11 */    04005,
	/* 12 */    010123,
	/* 13 */    020033,
	/* 14 */    042103,
	/* 15 */    0100003,
	/* 16 */    0210013,
	/* 17 */    0400011,
	/* 18 */    01000201,
	/* 19 */    02000047,
	/* 20 */    04000011,
	/* 21 */    010000005,
	/* 22 */    020000003,
	/* 23 */    040000041,
	/* 24 */    0100000207,
	/* 25 */    0200000011,
	/* 26 */    0400000107,
	/* 27 */    01000000047,
	/* 28 */    02000000011,
	/* 29 */    04000000005,
	/* 30 */    010040000007,
	/* 31 */    020000000011, 
	/* 32 */    00020000007 };  
	/* Really 40020000007, but we're omitting the high order bit */

/* Globals for powers of 2 for each
   value of w */
static int nw[33] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4), 
	(1 << 5), (1 << 6), (1 << 7), (1 << 8), (1 << 9), (1 << 10),
	(1 << 11), (1 << 12), (1 << 13), (1 << 14), (1 << 15), (1 << 16),
	(1 << 17), (1 << 18), (1 << 19), (1 << 20), (1 << 21), (1 << 22),
	(1 << 23), (1 << 24), (1 << 25), (1 << 26), (1 << 27), (1 << 28),
	(1 << 29), (1 << 30), (1 << 31), -1 };

/* Device and host globals -- 2^w -1 */
int *nwm1_d = NULL;

static int nwm1[33] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1 << 4)-1, 
	(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1, (1 << 9)-1, (1 << 10)-1,
	(1 << 11)-1, (1 << 12)-1, (1 << 13)-1, (1 << 14)-1, (1 << 15)-1, (1 << 16)-1,
	(1 << 17)-1, (1 << 18)-1, (1 << 19)-1, (1 << 20)-1, (1 << 21)-1, (1 << 22)-1,
	(1 << 23)-1, (1 << 24)-1, (1 << 25)-1, (1 << 26)-1, (1 << 27)-1, (1 << 28)-1,
	(1 << 29)-1, (1 << 30)-1, 0x7fffffff, 0xffffffff };

/* Kernel for multiplying into a given region, where original 
   region contains nothing. */
__global__ void mult_w8_null(unsigned char *ur1_d, 
		unsigned char *ur2_d, 
		int *multby_d, 
		int *prim_poly_d, 
		int *nwm1_d) { 
	// Calculate the index of the element.
	int index = blockIdx.x*blockDim.x + threadIdx.x;
	char x; 
	char y;
	int z;
	x = ur1_d[index];
	y = *multby_d;
	char prod;
	int i, j, ind;
	int k;
	char scratch[33];
	// Calculate the product from scratch
	prod = 0;
	for (i = 0; i < 8; i++) {
		scratch[i] = y;
		if (y & (1 << (7))) {
			y = y << 1;
			y = (y ^ prim_poly_d[8]) & nwm1_d[8];
		} else {
			y = y << 1;
		}
	}
	for (i = 0; i < 8; i++) {
		ind = (1 << i);
		if (ind & x) {
			j = 1;
			for (k = 0; k < 8; k++) {
				prod = prod ^ (j & scratch[i]);
				j = (j << 1);
			}
		}
	}
	// Place product result in region, where region 
	// where result is 0 if either original 
	// operand is 0.
	ur2_d[index] = (x == 0 || y == 0) ? 0 : prod;
}

/* Kernel for multiplying when the initial region is not null */
__global__ void mult_w8_notnull(unsigned char *ur1_d, 
		unsigned char *ur2_d, 
		int *multby_d, 
		int *prim_poly_d, 
		int *nwm1_d) {
	// Calculate the index of the element
	// to be operated on.
	int index = blockIdx.x*blockDim.x + threadIdx.x;
	char x; 
	char y;
	int z;
	x = ur1_d[index];
	y = *multby_d;
	char prod;
	int i, j, ind;
	int k;
	char scratch[33];

	// Code is exactly the same as previous kernel
	// except for last step.
	prod = 0;
	for (i = 0; i < 8; i++) {
		scratch[i] = y;
		if (y & (1 << (7))) {
			y = y << 1;
			y = (y ^ prim_poly_d[8]) & nwm1_d[8];
		} else {
			y = y << 1;
		}
	}
	for (i = 0; i < 8; i++) {
		ind = (1 << i);
		if (ind & x) {
			j = 1;
			for (k = 0; k < 8; k++) {
				prod = prod ^ (j & scratch[i]);
				j = (j << 1);
			}
		}
	}
	// The only difference between null and not_null
	// is that the product result is XORed with the
	// value currently in the region.
	ur2_d[index] = (x == 0 || y == 0) ? ur2_d[index]^0 : (ur2_d[index]) ^ prod;
}

/* Function called by benchmarks */
void galois_w08_region_multiply(char *region,      /* Region to multiply */
		int multby,       /* Number to multiply by */
		int nbytes,        /* Number of bytes in region */
		char *r2,          /* If r2 != NULL, products go here */
		int add)
{
	char scratch[33];
	unsigned char *ur1, *ur2, *cp;
	unsigned char prod;
	int i, srow, j;
	unsigned long l, *lp2;
	unsigned char *lp;
	int sol;
	cudaError_t ierr;
	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	// Check to make sure that primitive polynomial and
	// all other relevant data has been allocated
	// and copied onto GPU (device)
	if (prim_poly_d == NULL) {
		ierr = cudaMalloc((void **)&(prim_poly_d),sizeof(int)*33);
		cudaMemcpy(prim_poly_d, prim_poly, sizeof(int)*33, cudaMemcpyHostToDevice);
	}
	if (nwm1_d == NULL) {
		ierr = cudaMalloc((void **)&(nwm1_d),sizeof(int)*33);
		cudaMemcpy(nwm1_d, nwm1, sizeof(int)*33, cudaMemcpyHostToDevice);
	}
	unsigned char *ur1_d, *ur2_d;
	int * multby_d;

	// Pad the region to a multiple of 512
	// to take advantage of maximizing the number of threads.
	int total = nbytes+512-nbytes%512;

	// Allocate space for region if it has not been allocated
	// before and set globals.  Otherwise, grab the address
	// from the globals.
	if (c_ur1_d == NULL) {
		ierr = cudaMalloc((void **)&(ur1_d),sizeof(unsigned char)*(total));
		ierr = cudaMalloc((void **)&(ur2_d),sizeof(unsigned char)*(total));
		c_ur1_d = ur1_d;
		c_ur2_d = ur2_d;
	}
	else {
		ur1_d = c_ur1_d;
		ur2_d = c_ur2_d;
	}
	ierr = cudaMalloc((void **)&(multby_d),sizeof(int));
	// Copy regon and multiply_by number to device
	cudaMemcpy(ur1_d, ur1, nbytes, cudaMemcpyHostToDevice);
	cudaMemcpy(multby_d, &multby, sizeof(int), cudaMemcpyHostToDevice);

	// Call the kernels, error check, and copy results back 
	// onto host from device
	if (r2 == NULL || !add) {
		mult_w8_null<<<total/512, 512>>>(ur1_d, ur2_d, multby_d, prim_poly_d, nwm1_d);
		if (cudaSuccess != cudaGetLastError()) {
			fprintf(stderr, "Error!\n");
		}
		cudaMemcpy(ur2, ur2_d, nbytes, cudaMemcpyDeviceToHost);
	} else {
		cudaMemcpy(ur2_d, ur2, nbytes, cudaMemcpyHostToDevice);
		mult_w8_notnull<<<total/512,512>>>(ur1_d, ur2_d, multby_d, prim_poly_d, nwm1_d);
		if (cudaSuccess != cudaGetLastError()) {
			fprintf(stderr, "Error!\n");
		}
		cudaMemcpy(ur2, ur2_d, nbytes, cudaMemcpyDeviceToHost);
	}
	return;
}

/* Kernel for w = 16 and region is null */
__global__ void mult_w16_null(unsigned short *ur1_d, 
		unsigned short *ur2_d, 
		int *multby_d, 
		int *prim_poly_d, 
		int *nwm1_d) {
	// Calculate index of element to be operated on
	int index = blockIdx.x*blockDim.x + threadIdx.x;
	short x; 
	short y;
	int z;
	x = ur1_d[index];
	y = *multby_d;
	short prod;
	int i, j, ind;
	int k;
	short scratch[33];

	// Calculate product using primitive polynomial
	prod = 0;
	for (i = 0; i < 16; i++) {
		scratch[i] = y;
		if (y & (1 << (15))) {
			y = y << 1;
			y = (y ^ prim_poly_d[16]) & nwm1_d[16];
		} else {
			y = y << 1;
		}
	}
	for (i = 0; i < 16; i++) {
		ind = (1 << i);
		if (ind & x) {
			j = 1;
			for (k = 0; k < 16; k++) {
				prod = prod ^ (j & scratch[i]);
				j = (j << 1);
			}
		}
	}
	// Put the product in proper place in region
	ur2_d[index] = (x == 0 || y == 0) ? 0 : prod;
}

/* Kernel for w = 16 when the region contains data */
__global__ void mult_w16_notnull(unsigned short *ur1_d, 
		unsigned short *ur2_d, 
		int *multby_d, 
		int *prim_poly_d, 
		int *nwm1_d) {
	// Calculate index
	int index = blockIdx.x*blockDim.x + threadIdx.x;
	short x; 
	short y;
	int z;
	x = ur1_d[index];
	y = *multby_d;
	short prod;
	int i, j, ind;
	int k;
	short scratch[33];
	// Calculate product using primitive polynomial
	prod = 0;
	for (i = 0; i < 16; i++) {
		scratch[i] = y;
		if (y & (1 << (15))) {
			y = y << 1;
			y = (y ^ prim_poly_d[16]) & nwm1_d[16];
		} else {
			y = y << 1;
		}
	}
	for (i = 0; i < 16; i++) {
		ind = (1 << i);
		if (ind & x) {
			j = 1;
			for (k = 0; k < 16; k++) {
				prod = prod ^ (j & scratch[i]);
				j = (j << 1);
			}
		}
	}
	// XOR the product with the element and 
	// put the result in the proper place.
	ur2_d[index] = (x == 0 || y == 0) ? ur2_d[index]^0 : (ur2_d[index]) ^ prod;
}


/* Function called by benchmarks for w = 16 */
void galois_w16_region_multiply(char *region,      /* Region to multiply */
		int multby,       /* Number to multiply by */
		int nbytes,        /* Number of bytes in region */
		char *r2,          /* If r2 != NULL, products go here */
		int add)
{
	unsigned short *ur1, *ur2, *cp;
	int prod;
	int i, log1, j, log2;
	unsigned long l, *lp2, *lptop;
	unsigned short *lp;
	int sol;
	cudaError_t ierr;

	ur1 = (unsigned short *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned short *) r2;
	nbytes /= 2;
	
	// If primitive polynomial and other relevant
	// data are not already present on device,
	// allocate them and copy to device.
	if (prim_poly_d == NULL) {
		ierr = cudaMalloc((void **)&(prim_poly_d),sizeof(int)*33);
		cudaMemcpy(prim_poly_d, prim_poly, sizeof(int)*33, cudaMemcpyHostToDevice);
	}
	if (nwm1_d == NULL) {
		ierr = cudaMalloc((void **)&(nwm1_d),sizeof(int)*33);
		cudaMemcpy(nwm1_d, nwm1, sizeof(int)*33, cudaMemcpyHostToDevice);
	}
	unsigned short *ur1_d, *ur2_d;
	int *multby_d;
	int total;
	// Pad the region to a multiple of 512
	// to maximize the threads in use.
	total = nbytes+512-nbytes%512;
	// Allocate space for regions, if it is not already
	// allocated.
	if (s_ur1_d == NULL) {
		ierr = cudaMalloc((void **)&(ur1_d),sizeof(short)*total);
		ierr = cudaMalloc((void **)&(ur2_d),sizeof(short)*total);
		s_ur1_d = ur1_d;
		s_ur2_d = ur2_d;
	}
	else {
		ur1_d = s_ur1_d;
		ur2_d = s_ur2_d;
	}
	ierr = cudaMalloc((void **)&(multby_d),sizeof(int));
	// Copy region over to the device.
	cudaMemcpy(ur1_d, ur1, sizeof(short)*nbytes, cudaMemcpyHostToDevice);
	cudaMemcpy(multby_d, &multby, sizeof(int), cudaMemcpyHostToDevice);

	// Call kernels, check for errors, and copy the result back
	// to the host.
	if (r2 == NULL || !add) {
		mult_w16_null<<<total/512, 512>>>(ur1_d, ur2_d, multby_d, prim_poly_d, nwm1_d);
		if (cudaSuccess != cudaGetLastError()) {
			fprintf(stderr, "Error!\n");
		}
		cudaMemcpy(ur2, ur2_d, sizeof(short)*nbytes, cudaMemcpyDeviceToHost);
	} else {
		cudaMemcpy(ur2_d, ur2, sizeof(short)*nbytes, cudaMemcpyHostToDevice);
		mult_w16_notnull<<<total/512, 512>>>(ur1_d, ur2_d, multby_d, prim_poly_d, nwm1_d);
		if (cudaSuccess != cudaGetLastError()) {
			fprintf(stderr, "Error!\n");
		}
		cudaMemcpy(ur2, ur2_d, sizeof(short)*nbytes, cudaMemcpyDeviceToHost);
	}
	return;
}

/* Kernel for w = 32 when region does not have any data */
__global__ void mult_w32_null(unsigned int *ur1_d, 
		unsigned int *ur2_d, 
		int *multby_d, 
		int *prim_poly_d, 
		int *nwm1_d) {
	// Calculate index of element to be operated on.
	int index = blockIdx.x*blockDim.x + threadIdx.x;
	int x; 
	int y;
	int z;
	x = ur1_d[index];
	y = *multby_d;
	int prod;
	int i, j, ind;
	int k;
	int scratch[33];

	// Calculate product from primitive
	// polynomial.
	prod = 0;
	for (i = 0; i < 32; i++) {
		scratch[i] = y;
		if (y & (1 << (31))) {
			y = y << 1;
			y = (y ^ prim_poly_d[32]) & nwm1_d[32];
		} else {
			y = y << 1;
		}
	}
	for (i = 0; i < 32; i++) {
		ind = (1 << i);
		if (ind & x) {
			j = 1;
			for (k = 0; k < 32; k++) {
				prod = prod ^ (j & scratch[i]);
				j = (j << 1);
			}
		}
	}
	// Put product result in proper place.
	ur2_d[index] = (x == 0 || y == 0) ? 0 : prod;
}

/* Kernel for w = 32 when the region is non-empty */
__global__ void mult_w32_notnull(unsigned int *ur1_d, 
		unsigned int *ur2_d, 
		int *multby_d, 
		int *prim_poly_d, 
		int *nwm1_d) {
	// Calculate index
	int index = blockIdx.x*blockDim.x + threadIdx.x;
	int x; 
	int y;
	int z;
	x = ur1_d[index];
	y = *multby_d;
	int prod;
	int i, j, ind;
	int k;
	int scratch[33];

	// Calculate product using primitive polynomial
	prod = 0;
	for (i = 0; i < 32; i++) {
		scratch[i] = y;
		if (y & (1 << (31))) {
			y = y << 1;
			y = (y ^ prim_poly_d[32]) & nwm1_d[32];
		} else {
			y = y << 1;
		}
	}
	for (i = 0; i < 32; i++) {
		ind = (1 << i);
		if (ind & x) {
			j = 1;
			for (k = 0; k < 32; k++) {
				prod = prod ^ (j & scratch[i]);
				j = (j << 1);
			}
		}
	}
	// XOR product with original data,
	// and put result in proper place 
	ur2_d[index] = (x == 0 || y == 0) ? ur2_d[index]^0 : (ur2_d[index]) ^ prod;
}

/* Function called from benchmarks for w = 32 */
void galois_w32_region_multiply(char *region,      /* Region to multiply */
		int multby,       /* Number to multiply by */
		int nbytes,        /* Number of bytes in region */
		char *r2,          /* If r2 != NULL, products go here */
		int add)
{
	cudaError_t ierr;
	unsigned int *ur1, *ur2, *cp, *ur2top;
	unsigned long *lp2, *lptop;
	int i, j, a, b, accumulator, i8, j8, k;
	int acache[4];

	ur1 = (unsigned int *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned int *) r2;
	nbytes /= sizeof(int);

	// Allocate primitive polynomial and other relevant 
	// data if they have not already been allocated and
	// copied onto the GPU (device).
	if (prim_poly_d == NULL) {
		ierr = cudaMalloc((void **)&(prim_poly_d),sizeof(int)*33);
		cudaMemcpy(prim_poly_d, prim_poly, sizeof(int)*33, cudaMemcpyHostToDevice);
	}
	if (nwm1_d == NULL) {
		ierr = cudaMalloc((void **)&(nwm1_d),sizeof(int)*33);
		cudaMemcpy(nwm1_d, nwm1, sizeof(int)*33, cudaMemcpyHostToDevice);
	}
	unsigned int *ur1_d, *ur2_d;
	int *multby_d;
	int total;
	// Pad region size to a multiple of 512
	// to maximize thread use.
	total = nbytes+512-(nbytes)%512;
	if (l_ur1_d == NULL) {
		ierr = cudaMalloc((void **)&(ur1_d),sizeof(int)*total);
		ierr = cudaMalloc((void **)&(ur2_d),sizeof(int)*total);
		l_ur1_d = ur1_d;
		l_ur2_d = ur2_d;
	}
	else {
		ur1_d = l_ur1_d;
		ur2_d = l_ur2_d;
	}

	ierr = cudaMalloc((void **)&(multby_d),sizeof(int));
	// Copy region onto GPU
	cudaMemcpy(ur1_d, ur1, sizeof(int)*(nbytes), cudaMemcpyHostToDevice);
	cudaMemcpy(multby_d, &multby, sizeof(int), cudaMemcpyHostToDevice);

	// Call kernel, error check, and copy result from device to host.
	if (r2 == NULL || !add) {
		mult_w32_null<<<total/512,512>>>(ur1_d, ur2_d, multby_d, prim_poly_d, nwm1_d);
		if (cudaSuccess != cudaGetLastError()) {
			fprintf(stderr, "Error!\n");
		}
	} else {
		cudaMemcpy(ur2_d, ur2, sizeof(int)*(nbytes), cudaMemcpyHostToDevice);
		mult_w32_notnull<<<total/512,512>>>(ur1_d, ur2_d, multby_d, prim_poly_d, nwm1_d);
		if (cudaSuccess != cudaGetLastError()) {
			fprintf(stderr, "Error!\n");
		}
	}
	cudaMemcpy(ur2, ur2_d, sizeof(int)*(nbytes), cudaMemcpyDeviceToHost);
	return;
}
\end{verbatim}
\end{tiny}
